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EXTRAPOLATION TECHNIQUES APPLIED TO MATRIX METHODS IN 
NEUTRON DIFFUSION PROBLEMS ^ 


By Robert R. McCbbabt 


SUMMARY 


SYMBOLS 


A general mairix method is developed for the ' solution of 
characteristic-value problems of the type arising in many 
physical applications. The scheme employed is essentially that 
of Qauss and Seidel with appropriate modifications needed to 
make it applicable to characteristic-lvalue problems. An itera- 
tive procedure produces a sequence of estimates to the answer; 
and extrapolation techniques, based upon previous behavior of 
iterants, are utilized in speeding convergence. Theoretically 
sound limits are placed on the magnitude of the extrapolation 
that may be tolerated. 

This mairix method is applied to the problem of finding 
criticality and neutron fluxes in a nuclear reactor with control 
rods. The two-dimensional finite-difference approximation to 
the two-group neutron-diffusion equations is treated. Results 
for this example are indicated. 

The calculations were peiform^ on the IBM card-programmed 
calculator. 


INTRODUCTION 

A general matrix metbod is developed for the solution of 
characteristic-value problems of a type arising in many 
physical applications. The method of this paper is essen- 
tially tliat of Gauss and Seidel (ref. 1), which itself is but a 
special case of the method of conjugate gradients (ref. 2). 
The adaptation of the Gauss-Seidel technique to the charac- 
teristic-value problem calls for a means of computing suc- 
cessive estimates of the characteristic value as weU as the 
vector. This calculation is made to rely upon Tmrner’s 
technique (ref. 3) for assigning a meaning to the ratio of two 
vectors. 

Extrapolation techniques are also employed to speed up 
the convergence of the iterative process. One of these is 
based on Turner’s original formula (ref. 3), and the other is 
a slightly more complicated modification. 

The number of iterations required for convergence is not 
studied theoretically here as in the “a-step” methods, but 
the minimization of a suitable form at each step is derived. 

The method is applied to two-group neutron-diffusion 
equations. The calcidations were performed at the NACA 
Lorvis laboratory. 


The following symbols are used in this report: 


A, B, L, U 
Bl 

D, E, F, G, J, X 
h 

i, j, k 

koi 

N 

Pv, 

r 

Tc 

sgn a= 
t 

W,a> 

T 

At‘^1 

X 

T 


matrices 
axial leakage 
vectors 

grid dimension 
indices 

thermal multiphcation constant 
average square slowing down length for 
fast neutrons 

average square diffusion length for ther- 
mal neutrons 
number of nuclei per cc 
resonance escape probability 
radial coordinate 
core radius 
[ 1 a>0 
-1 a<0 
1 0-a=0 

reflector thickness 
weight functions 
characteristic value 

deviation at ith point of /rth iteration (eq. 
(80)) 

difference X'** 

actual damping rate 
bulk damping rate 
neutron fluxes 


Parameter groupings: 


a=~+Bi 


1 

a j2 


b= 


'tTt jJlo 




c=7^+-B| 

„ 1 

^tT, tJtQ -LifQ 


/=4;+B! 


J-'tbn 




• Supersedes NAOA TN 3611, “Extrapolation Techniques Applied to Matrix Methods In Nentron DJlIasion Problems," by Robert R. McCready, 1966, 
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Subscripts: 

f 

th 

tr 

0 

1 

2 


■^thi 


m= 




\ Pt^i T 1 

^tr, l4j 


fast 

thermal 

transport 

reactor 

reflector 

rod 


■L'th, 


THE METHOD 

MATBIX FORMULATION 

Consider the matrix equation 

AX=yBX (1) 


where A and B are nY.n matrices, X is an ii-component 
vector, and the characteristic value of 7 is a scalar to be 
determined. A may be separated into the sum of two 
triangular matrices L and U, where L contains all the 
diagonal elements of the original matrix A. 

This separation, which anticipates the Gauss-Seidel 


process, is effected in the following manner: 

A=L+U ( 2 ) 

j<i; lti=0 j>i (3) 

«w=ao j>i; ^,{,=0 j<i (4) 

If i is a nonsingular matrix (always true if for all i), 
equation ( 1 ), modified to 

(L+U)X=yBX (5) 

may be multiplied bj^ giving 

{I+L-^U)X=yL-^BX ( 6 ) 


For a given X, the quantities L~^UX and L~^BX of 
equation ( 6 ) may be calculated without the actual formation 
of L~^. This fact, which is very helpful for systems con- 
taining matrices, arises in the following manner and de- 
pends upon the triangular nature of L. Let D be the 
vector defined by 

D=L-^UX (7) 

Then 

LD==UX=C ( 8 ) 

whence 

^iidi=Ci (9) 

gives di, since aU the c* can be computed from U and X in 

equation ( 8 ). Then 

f2ldl~t-Z22dj = Ca (10) 

gives dj, and 

hidiAhiditAlsida^Cz ( 11 ) 


ITERATIVE SCHEME 

Equation ( 6 ) may be written 

X=yL-^BX-L-^UX ( 12 ) 


which may be inteipreted as defining the iterative scheme 


X,+,=y,+rL-^BX,~L-^UX, (13) 

in which 7 t +2 is an estimate to 7 that can be calculated from 
Xf To obtain 7 *+i, form the inner product of the vector 
sgn L~^BXt with each side of equation ( 6 ) ; thus. 


.. (sgn L-^BX„(I+L-^U)X,) 
(sgn L->5Xt,L-‘BX,) 


(14) 


Equations (13) and (14) are the basic equations of the 
iterative scheme. Given any Xt, 7 *+i is computed from 
equation (14) and 7^+1 and Xt are placed in (13) to 3 dold 
the next iterant This process is repeated until X* 

and 7 t+i converge. 

Some normalization is necessary in problems of a homo- 
geneous nature. The simplest method of normalization is 
to adjnst a permanently specified coordinate of X* to unity 
before beginning each iteration. This is accomplished 
by dividing each element of the vector by the specified 
coordinate. 

The ratio defined by equation (14) was chosen for sim- 
plicity of calciilation on available punched-coi’d equipment. 
That ratio can be compared to the Kayleigh quotient 
(for eq. (13)) 


y _(AA) 


(16) 


where 


Jt=L-^BXt 


(16) 


0,= (I+L->U)X, (17) 


by noting that each of the relations (14) and (16) constitutes 
a weighted smn of local (point by point) values y^^.\ of yi+i. 
These local values are defined by 

(18) 

where and are the zth components of Ot and J*, 
respectively. The weighted average associated Avith (16) is 

7t+i=S'a^(i't+i (19) 

i 

where 

""'""soft 

m 


while the weighted average associated with (14) is 

i 

where 


£0,= 


\i 




( 21 ) 

( 22 ) 


gives dj, and so forth, so that L~^ need not be computed in 
order to obtain L~^UX. The same ai-gument applies 
to L-^BX. 


Equation (15) selects that value of 7^+1 which minimizes 
the sum of the squai’es of the residuals of equation ( 6 ) when 
that quantity is thought of as a function of 7 ^+ 1 . The sum 
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of tli 6 squares of the residuals is not, of course, the only 
quadratic form that is suitable for minimization (ref. 2 ). 
Consider the expression 


wliich is zero for X^=X and 7 t+i= 7 . This generally positive 
quantity can be minimized (made closer to its ultimate value 
zero) by sotting 


(sgn 

(sgnJt,Jt) 


(24) 


which is equation (14) expressed in terms of J and 0. 


EXAMPLE 

To illustrate the convergence of this method in a special 
case, consider the problem of equation ( 1 ) with 



evaluate a bulk damping rate that describes in an average 
way the over-all trend of the individual components of the 
iterant vectors. 

Assume that each iterant X^ is made up of the sum of the 
solution X and two error vectors Et and Ft satisfying the 
damping relations 


and 


One may compute 


Et+i-^-rEt 

(27) 

Ft+i= — rFt 

(28) 

relations hold : 


Xo=X-\-Eo-\-Fo 

(29) 

Xi—X-{-tEo — tFo 

(30) 

X3=X-\-FEq-^FFo 

(31) 

A 3 = A-k FEo - FFo 

(32) 

, A 3 Aj 
^~Ai-Ao 

(33) 


The “vector division” implied in equation (33) is possible 
because, under the initial assumption of error behavior 
(eqs. (27) and (28)), the vectors X^—Xi and Xo are 
collineax and therefore differ only in length. 

If the error vectors are eliminated from equations (30) 
and (32), one obtains 

( 34 ) 


which has the real solution = 1.020070, = 1.329658, 

1.000000, = 1.1 09886; 7=0.549429 and two solu- 

tions wth complex characteristic values. This solution was 
foimd by the ordinary process of solving the characteristic 
equation. 

This problem was solved in 15 iterations starting with an 
initial guess of Ao= (10,100,1,1000). The values of suc- 
cessive iterants, together with those of 7 , are listed in the 
following table. The iterants are normalized so that 
AT’ = 1 at the start of each iteration : 


which gives the answer as a linear combination of the alter- 
nate iterants Ai and A 3 . 

The preceding analysis suggests that a formula analogous 
to (34) be used to estimate the answer. The difficulty here 
is that equation (33) may be meaningless when equations 
(29) to (32) do not hold. To circumvent this difficulty, a 
method of computing F is needed. Toward this end, 
define 5t+i by means of 

(35) 


k 

x<» 


x<« 

k 

V. 

0 

10 

100 

1000 


1 

-a 692996 

1. 069635 

L 023212 

-la 634689 

2 

.271412 

.438316 

.812771 

-1.362363 

3 

.848870 

L 666162 

1.188384 

-.097105 

4 

U 132470 

1.360721 

L 120240 

.626664 

6 

L 006120 

L 292707 

L 097668 

.892807 

6 

L 010639 

1.337044 

1. 112348 

.638366 

7 

1.024644 

1. 331918 

LU0639 

.647271 

8 

L 019872 

L 328104 

L 109368 

.661383 

0 

1.019603 

L326848 

L 109049 

.649107 

10 

1.020236 

L329788 

1. 109929 

.649291 

11 

1. 020078 

1.329696 

1.109866 

.649610 

12 

1. 020048 

1.329660 

L 109886 

.649423 

13 

1. 020076 

1.329666 

1.109889 

.649422 

14 

1. 020071 

1.329666 

1.109886 

.649432 

16 

L020070 

1.329668 

1.109886 

.649430 


EXTEAPOLATION TECETNIQUE 

If, instead of four components, the iterant vector has many 
components, techniques of extrapolation are usually desirable 
to speed convergence of the process. The technique em- , 
ployed here, which is due to Turner (ref. 3), attempts to 


and define F by means of 

2D5T sgn 
J 

The direct analogy to equation (14) will be noticed. Equa- 
tion (36) permits computation, in on average way, of the 
damping of the error vectors. With F available, the extra- 
polated value A' of A is computed from 

(37) 

In case the error is damping exactly as assumed in (27) and 
(28), equation (36) gives the value indicated by (33), and 
equation (37) reduces to (34); that is, X' becomes the 
answer A. 

’ Since the ideal damping behavior is rarely an actuality, it 
is of interest to examine the effect of the preceding process on 
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error components. Suppose that Xo is more adequately 
represented by 

Xo=X+±Ei^ (38) 

f=i 


S(X,t)= 


" 1-T* 


J 

1-^T» 


(47) 


where has a damping rate (positive or negative) of Xj. 
Then 

X,=X+±.\}Ei*^ (39) 

and 

(40) 

j-i 

hold. The extrapolation indicated in equation (37) now 
yields the following relation between the estimate X' and 
tlie answer X: 

X'=X+J: Ei‘> (41) 

(-1 1 — r 

This interpretation is useful, since it indicates the damping 
effect upon the errors of three iterations and one extrapo- 
lation. 

If, for simplicity, one of the errors E^^ and its damping 
rate X| are designated by E and X, respectively, then 


In formula (47), j>4=. The second factor places a zero 
(maximum damping) at just the points of minimiun damping, 
that is, at the values of X determined by (44). If now 
dB/d\ is taken as zero and the limit ± 1 is placed upon the 
resulting RatiT^), the limiting safe values of t* are obtained 
by finding the least of the roots of 

[jr^—2{j—iy-‘^->r{j—2y-*\r^±[j—2{j—\)T^-\-{j—4.)T*\=Q 

(48) 

where r* satisfies 

j., (i-l)(i-2)dby5/-12i-H4 

The revised formula (47) has both the effect of ensming that 
no component will be impaired in its damping by the extra- 
polation and also that the least rapidly damping component 
receives a zero contribution in the extrapolation. 

Since, as before, for some error component E, 

Xi=X+\^E (60) 


i?(X,r) 


(X*-r*)X^-’ 

l-r" 


(42) 


gives the damping of this error component as a result of j 
iterations and one extrapolation. The “extreme” value of 
R (actually that value farthest from zero; i. e., farthest from 
maximum damping) may be foxmd by setting 


This yields 


di? 

1 

r 

■r^(i-2)V-3 „ 

(43) 

dX" 


1-r^ " 



'•7— 2X 



x^=( 

. J J 

(44) 


as the equation to be solved for the values of X which are asso- 
ciated with the errors that receive the minimum damping 
from the process of j iterations and one extrapolation. 
Equations (42) and (44) give Rai{j), the extreme value of 
jB, as a fimction only of t and j: 



To find the value of -ri so that this function (K«i,) cannot 
exceed the bounds ±1 [i. e., so that the slowest damping 
component (and hence all components) cannot increase 
through extrapolation], r must be less in absolute value than 
the least of the roots of 

If only such are used, the convergence of the process 
cannot be impaired by the extrapolation. 

Suppose now that the previous value of R (X,r) is replaced 
by the formula 


X,_a=X-|-X^-*E (61) 

Xi.i=X+V-*E (62) 


in which X represents the answer, the specification of (47) 
as a damping formula imphes 


^ , jV-2(i-l)TV-»-Ki-2)TV-^ „ 
j-2{j-\y+{j-2y 


(63) 


where X' is the extrapolated value of Xj. If >JE, }d~^E, and 
\^~*E are eliminated from (63) using the relations (60), 
(51), and (52), then 


jx,-2{j-iyx,_,+{j-‘^)r*x,_, 

j-2{j-iy+{j-2y 


(64) 


Comparison of (42) (with j=3) with (37) on one hand and 
of (47) with (54) on the other hand leads to the following 
valid rule of thumb to obtain the extrapolated value of X 
for a given damping fimction: Replace the power X' of X 
in the damping function by Xi ; the resulting linear combina- 
tion of alternate iterants is the formula for the extrapolated 
X. It is easily verified that the validity of this arises from 
the manner in which the error vectors are assmned to behave. 

The smallest roots of equations (46) and (48) are listed 
in the following tables: 


Eq. (48) 

S 


4 

0.6667 

6 

.8746 

8 

.0213 

10 

.W27 


Eq. (46) 

} 

T* 

4 

aS284 

6 

.8011 

8 

.0233 

10 

.0300 


These are the upper limits of the “safe” values of within 
the framework of the definition. 
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APPLICATION TO REACTOR THEORY 

GENBRAI. REMARKS 

Multigroup reactor equations can be solved, in principle, 
by the present method. The number of components in the 
vector solution, to be discussed in detail later, is approxi- 
mately equal to the product of the number of grid points 
and the number of groups in the multigroup scheme. An 
extreme increase in the number of these elements lengthens 
the problem considerably. The calculations hero are per- 
formed in accordance with two-group neutron-diflFusion 
theory. 

The two-dimensional reactor with control rods, which is 
considered later, is suited to two-group calculations, since 
the control rods are particularly effective on the thermal 
group, and two-group calculations are good for thermal 
assembhes. 

The following illustration is introduced to show the general 
principles of the matrix setup in detail. These principles do 
not change for the more complicated two-dimensional prob- 
lem that is treated later. A relatively simple one-dimensional 
problem has been chosen to illustrate the detailed setup and 
the consequent matrix. 

EXAMPLE OF TWO-GROUP DIFFD8ION EQUATIONS 

The one-dimensional diffusion equations for a reflected 
tliormal reactor of slab geometry ore (ref. 4) 

(55) 

and 

cV(»+d¥v=0 (56) 

for the core, and 

^'-/V=0 (67) 

and 

Vth+'m<pf=Q (58) 

for the reflector. 

The parameters a', b, c', d, gr', and m are defined in the 
list of symbols; y is the characteristic value of the system and 
equals 1 for criticality. When y converges to a value other 
than unity, the uranium concentration is adjusted and the 
proeess repeated. 

The differential equations (55) to (58) are replaced by 
finite-difference equations; the operation dV/di* is estimated 
by means of the approximate formifla 

(59) 

where the points of the region are niunbered in order as grid 
points of a finite-difference net, and h is the distance between 
successive points. In the following, is the core radius, 
Vc+t the complete reactor radius, and point 6 lies on the 
interface: 


. 1 . 


0123456789 
0 r, 

The boundary conditions are that the fast and thermal 
fluxes have zero current across the plane of symmetry (a;=0). 
This condition 


=0 

dx 


can be approximated by 


for both ipf and tpa,. The condition of continuity of currents 
at the interface is met by approximating the derivatives in 
the expression 

-Vo^(r-0) = -X«,j^(r,-l-0) (61) 

for both the fast and thermal fluxes. The remaining condi- 
tion is that the flux be zero at the outer boimdary. If the 
fast flux is designated by <p and the thermal by rp, the system, 
becomes: 

Equation 






<Pi Vs X Vi Vs 

T,r./0 ^ ^ 


y8+ya~2v>7 


—fvi=^ 


g |/'9:s=0whichincorporates959=0 

n 1 


for the fast-balance equations and, for the thermals: 
Equation 


^1— ^0=0 


'!i±!±ip!=^-c%+dt,=o 

. 'I'B—'Ps ^ i'T—'f'S 


7 

8 which incorporates ^.=0 
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The variables ^ to <p 8 and to ^3 may now be written as Xo to 
Xi and X 9 to A'lr, respectively. The matrix formulation of 
these equations is presented in figure 1 . The following 


symbols have been introduced: 

(64) 

^=2/?t»+a (65) 

Fn=2!h?+f ( 66 ) 

Cn=2!h?+c (67) 

Q,=2l¥+g ( 68 ) 

L=^ (69) 

\r,fl 

(70) 

“ir. tAl 



FiatTBB 1. — Matruc formulation of equations (62) and (63). 


EECAPITOLATION 

To review the general appbcation of the method to two- 
group reactor equations, consider the following broad out- 
line of this process: 

(1) Write the two-group equations with the parameter T 
introduced as a multiplier of the production term of the fast- 
balance equations. 

( 2 ) Express the difl^erential equations by their finite- 
difference approximations so that they become a linear 
algebraic set of the type associated with equation ( 1 ). 

(3) Perform such iterations and extrapolations as neces- 
sary to obtain weU-converged values of T and X. 

(4) Adjust the rnmnum concentration and repeat step (3) 
using the original answer from (3) for the initial guess Xo. 
The concentration should be changed so that 7^1. 

( 5 ) Kepeat (3) and (4) until 7 converges. If criticality is 
desired, change the concentration so that the conveiged 
values of 7— >1. 


TWO-DIMENSIONAL REACTOR WITH CONTROL RODS 

GEOMETBY OF HEACTOE 

The reactor (see fig. 2 ) is cylindrical and water-reflected 
with a core composed of aluminum, water, and m-anhim, 
which are assumed to be homogeneously mixed. The 
height of the reactor is 70 centimeters, the outside radius 50 
centimeters, and the core radius 32 centimeteis. Nine 


cadmium control rods are inserted in the core; one, a cjdin- 
drical rod of 2 -centimeter radius, is centered along the axis 
of the reactor. The remaining eight rods are equally spaced 
on a radius of 24 centimeters and are shaped so as to be 
bounded by coordinate surfaces. Each of these rods extends 
over a radial distance of 4 centimeters and subtends a central 
angle of 9°. 

The symmetry of this assembly is an important factor in 
making solution of the reactor problem feasible. The flux 
in the 45° sector indicated in figme 2 is adequate to represent 
the flux in the entire reactor; in fact, additional symmetry 
within the sector implies that only half the sector need be 
considered. The three-dimensional problem is made two- 
dimensional (computation-wise) by estimating the neutron 
leakage in the axial direction due to the finite height of the 
reactor. This is based upon an axial cosine distribution 
similar to the bare pile solution (eq. (75)). 

COMPOSITION AND NUCLEAE PAEAMETEES 

The core volmne is proportioned between the water 
(density, 1 g/cc) and aluminum by assuming a volume ratio 
of aluminum to water of 0.76. The nuclear diffusion con- 
stants for the core and reflector are listed in the following 
tables. The subscripts 0 , 1 , and 2 refer to the core, reflector, 
and rod regions of the reactor, respectively: 



Parameters for the rod regions are imnecessary because of 
the simplified treatment of the rod, in which the thermal 
neutron flux is assumed to vanish on the rod boundary and 
the radial and axial leakages are assumed to balance in the 
absence of fast-neutron absorption processes. The thermal 
parameters in the preceding table are those associated ivith 
an atom ratio of N^/N" of 350; these, of course, change for 
different uranium concenti’ations. 



Figure 2. — Two-dimensional reactor. 
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EQUATIONS AND BOUNDARY CONDITIONS 

The two-group equations (ref. 4) for the core are taken to 
be 

VV/o — (^7=2 — ^/t)= — 7 ^(*0 yj~ (71) 

\L/0 / ■L'tho 

and 

VVrto — (yi — \~Bl \ Pm yr- <P/n=0 (72) 

V-O/Ao / ao L//0 

All the parameters of equations (71) and (72) ai-e the ordinary 
nuclear ones, except the ai-bitrarily inserted y, which is a 
measure of the criticality and is equal to 1 for a critical 
assembly. 

In the reflector the two-group diffusion equations take the 
form 

vV/i-(:^+5i) <fin=0 (73) 

VV«l — (yi — \ ytAi PiM yy <Pj\=(i (74) 

while the fast-diffusion equation for the rods is taken to be 

vVr2~-Bfv’/i=0 (75) 

Any change in rod boundary conditions would not affect the 
general piinciples of the numerical scheme. As will be seen 
from the boundary conditions, the thermal neutrons do not 
require a diffusion equation within the rods. The region 
considered in the problem is a 22%° sector (of the circle of 
fig. 2), one side of which passes through the center of one of 
the outlying rods. This is illustrated in figure 3. The 
symmetiy of the over-all reactor imphes that the normal 
derivatives of the flux across the surfaces A and B are zero. 
This implies that the fiux at aU points of the circle of figure 2 
can be foimd by solution for the fiux only in the sector 
indicated in figure 3. The condition of continuity of fluxes 
and currents is involved at the core-reflector interface, in- 
dicated by O' in figm-e 3. The vanishing of the fast and 
theimal flux at the outer boundary (Z?) is also required. The 
thermal flux is taken as zero on the edge of the control rod, 
and the continuity of the fast flux and cmTent is considered 
to hold on the core-rod interfaces. The details of the 
mathematical fonnulation of these conditions are deferred 
imtil the general discussion of the difference equatiops. 

FINITE-DIFFERENCE EQUATIONS 

In order to write the reactor equations as finite-difference 
approximations, the sector of figm’e 3 is divided into a grid 
net of points. The flux is determined by solution of the 


Radius, 

cm 



Figure 4. — Reactor grid points. 


linear algebraic system of equations that results from writing 
the finite-difference approximation to the fast- and thermal- 
diffusion equations at each point. The grid arrangement 
used in the present problems is indicated in figure 4. The 
thermal flux (components 1 to 139 of the vector solution) has 
the following breakdown into groups of components: reflector 
(1 to 73), core-reflector interface (74 to 79), core (80 to 
139). The fast flux (140 to 291) has the following break- 
down: reflector (140 to 212), core-reflector interface (213 to 
218), core (219 to 284), control rods (231, 232, 237, 238, 243, 
244, and 285 to 291). 

The number of components associated with the thermal 
and fast fluxes, respectively, differs because of the assign- 
ment of boundary conditions at the control rods, which 
brings the fast flux into a larger area of definition. 
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For two-diinensional cylindrical geometry, tlie operation 


of the form vV is given by 


dr ■’■ 5 ^ 


(76) 


This form is to be replaced by a difference opeitition that 
relates each point to its four nearest neighbors. If the point 
in question is designated by the subscript zero and the others 
are 

1 


hr 

4: • ha 0- he ■ 2 r | 6—* 
hr 

3 


where hr and he are the grid widths in the r and 6 directions, 
respectively, then at <p=<fXi the following approximation is 
used: 





(77) 


With this designation (and barring certain exceptional 
points to be discussed later), one may move from point to 
point on the grid and write equations of neutron balance 
for each of the two neutron groups. 

The foUoiving equations may be taken as typical illustra- 
tions: 

Thermal-balance equation 93 (see fig. 4): 


(e+2s) - 




2,2 1 


(78) 

/ A/r. «A0 -^/O 


Fast-balance equation 234: 




In contradistinction to equations (78) and (79), there are 
certain special equations that hold at the exceptional points 
referred to earlier. These equations result from one or more 
of die following conditions: 

(1) Continuity of currents at interfaces 

(2) Zero flux at the outside boimdary 

(3) Zero current across planes of symmetry 

(4) Change in grid dimensions 

Condition (1) is tirnted by matching a suitable ratio of 
normal derivatives from either side of the interface. Each 
of these derivatives is evaluated by a five-point differentia- 
tion formula. Condition (2) is treated by writing the 
difference approximation to the diffusion equation for points 
adjacent to the outside boundary with zero replacing the 
flux at the boundary point. 


Condition (3) is accmmted for by A\oiting the diffusioi 
equation of a point on the plane assuming the same flin 
at grid points on either side of the plane. For condition (4) 
0-wise interpolation formulas are used to define fluxe. 
at the points marked X on figure 4, and these are utilized 
where needed, in writing diffusion equations in the finer net 
If each equation of the set is written in order and the pro 
duction terms are isolated as illusti-ated in equation (71) 
then the rnatrix equation constructed from the approximat( 
fin ite difference may be written in the form of equation (1). 

. The matrix B is singular, largely consisting of zero ele- 
ments with an essentially diagonal group of nonzero temu 
somewhat off the leading diagonal. The matrix A has i 
substantial number of nonzero elements crowded quite 
dose to the leading diagonal. This latter situation is nu 
merically desirable, as elements far from the leading diagona 
tend to slow the convei^ence of numerical processes. 

If criticality is desired, the concentration of fissionable 
material is adjusted, after 7 and X have convei-ged, and a 
whole new set of calculations is nm until a new value foi 
7 is reached. This process may be continued until 7=1. 

The method can also be used to compute reactivity 
changes; the calculation time is again shortened consider- 
ably if flux distributions are not demanded. 

SOLUTIONS OF TWO-DIMBNSIONAL PROBLEM 

The results of the calculations of the supercritical 
(7=0.948) case are shown in figures 5 to 7. Figure 5 give, 
the fast flux as a function of r for 0=0°, 9°, and 18°. The 
control rods have no substantial effect on the fast flux. 
Figure 6 gives the corresponding thermal flux and shows the 
localized effect of the control rods. Figure 7 presents iso- 
flux contours of the thermal flux. The 0.19 contour in the 
reflector and the 0.234 contour in the core represent relative 
maximums. 

COMMENTS ON APPLICATION OP THE METHOD 

A number of numerical quantities may be examined in an 
attempt to evaluate the degree of convergence of a system. 
One of the most natural of these quantities is the sum of the 
squares of the residuals. Another may be formed by con- 
sidering the fact that, as the limit is approached, the ratio 
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Fiouhb 6. — Thermal-neutron flux for azimuth angles of 0°, 9°, and 18°- 

Tt+i/Tt+i must tend toward unity. This means that the 
deviation defined by 

(80) 

~t+i 


must tend toward zero. The average absolute value of the 
deviation, summed over all points of the reactor, is 


it+il 


n 

s 

1 

• 

yt+1 

n 


(81) 


where n is the number of reactor points. 

An illustration of the behavior of this quantity is given in 
the following table: 


i+l 

ua.i 

ha.u.. 

6 

0.02199 

0.2133 

IB 

.02661 

.3156 

23 

.002M 

.0147 

31 

.00130 

.ous 

40 

.000793 

.0061 

49 

.00028 

.0027 


The iterants listed are those which just precede the extra- 
polation process. These are chosen so as to minimize the 
effect of fluctuations introduced by the extrapolation 
technique. 

These illustrative values come from the second general 
process; that is, after y had converged to 1.2064, the con- 
centration (and hence elements of the matrices A and B) 
was changed and a new series of iterations begim. This 
converged (more rapidly than the first run) to a value of 0.948. 

^To estimate the value of uranium concentration needed for 
the new run, the equation 

1.2064 ka, (old)=^tt (new) (82) 

was used to compute a new koi from which to obtain a new 
concentration. This formula is an approximation, since the 
influence of a change in concentration upon if* is appre- 
ciable. The better rate of convergence of tbe second run is 
duo to the fact that the flux is relatively independent of the 
characteristic value, so that the initial estimate for the second 
run was a relatively good one. 

The quantity |Ai(}.i| reflects the convergence of y, which is 
faster than that of the vector X. 

In order to determine the degree of convei-gence of X, 
consider the quantities 


dt+i=g 1 (83) 

(84) 

and the maxim u m designated by Typical 

values of these quantities are as follows; 


y~vm 

k+l 


lai 


lU 

aoo90 

a000034 

0.000269 

122 

.0083 

.000028 

.000176 

133 

.0010 

.000017 

.000167 

144 

.0028 

.000010 

.0000601 

155 

.0023 

.000008 

.0000709 


T«0.948 

t+i 

<U+i 



6 

0.3209 

0.001137 

a 00918 

16 

.2139 

.000735 

.00530 

22 

.0223 

.000077 

.000201 

31 

.0044 

.000016 

.0000666 

40 

.0020 

.000007 

.0000216 


The sums of the squares of the residuals for the two cases 
■y= 1.206 and 7=0.948 are as follows: 


'1 

=IJ206 

k+l 


111 

9.44XlCr7 

122 

4.98X10-? 

133 

2. 66X10-? 

144 

8. 60X10-* 

155 

6.37XIO-* 


1 

=0.948 

t-R 

27?J 

6 

L63X10-* 

15 

3.51X10^ 

22 

1.80X10-* 

31 

1.22X10-' 

40 

L99X1(H 



Figtob 7. — C!ontour lines for tliermal flux. 
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Several general observations can be made about tbe process: 

(1) The number of iterations in this problem starting from 
an initial guess to a weU-converged value of X was about 
150 to 175. 

(2) In general, 8 to 10 iterations between extrapolations 
seem desirable, as the use of too few iterations does not allow 
the establishment of a fairly uniform damping rate. 

(3) The extrapolation formtda of equation (37) seems best 
for rough estimates where error components are being damped 
rapidly; that of equation (54) seems to be superior for latex- 
extrapolations where one is closer to the solution. 

(4) TVhen computed values of exceed the upper limit, 
they may be replaced by the limit from the tables giving • 
roots of equations (46) and (48) and then the extrapolation 
may be carried out, or two more iterations performed with 
T* recomputed xmtU it falls within presciibed limits. 

The following table gives the sum of the squares of the 
residuals for (a) direct iteration from Xn to (b) eight 
iterations from Xji followed by extrapolation with ‘V safe” 
when "t* computed” was too large, then iteration to 
(c) eight iterations from Xn followed by two iterations and 


a test imtil ‘V computed” was less than ‘V safe,” tlxon 
extrapolation followed by iteration to Xt^: 


Case 

2r; 

1 

1.3SX10-' 

3.44X10-* 

1.6SXKH 


Lewis Flight Propulsion Laboratory 

National Advisory Committee for Aeronautics 
Clevel.and, Ohio, May IS, 1956 
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